*Figure 2

cd ""

use "master_data_adults_child_s13456_shared.dta", clear 

drop meal_control snack_control
gen meal_control=1 if treatment==2
replace meal_control=0 if treatment==1
tab meal_control

gen snack_control=1 if treatment==3
replace snack_control=0 if treatment==1
tab snack_control

drop if type_id=="X" // **Dropping newborns born during the three years of study
drop if id=="3885553B" //**error in the id number and empty data
drop if id=="3885553G" //**error in the id number and empty data
replace treatment=2 if id=="38855553B" //**retrieving the correct treatment for this id
replace treatment=2 if id=="38855553G" //**retrieving the correct treatment for this id
drop if id=="101aa22c" //**this id does not exist
drop if id=="27655663" //**outlier on many measurements

***percentile in children
gen pzheight=normal(zht)
gen pzweight=normal(zwt)

*BMI
*Calculate z-scores for BMI
*Calculate z-scores for BMI
egen gendermax=max(gender), by(id_nb)
ge sex=(gendermax==1) 
replace sex=. if gendermax==.

drop zbmi pzbmi
egen zbmi = zanthro(bmi,ba,UK), xvar(months) gender(sex) gencode(male=1, female=0) ageunit(month)
gen pzbmi=normal(zbmi)

********************************************************************************
 ***Figure 2: Changes in the distribution of the BMI percentile rank in childre****
 
 //ksmirnov test//

gen session1_3=1 if session==3
replace session1_3=0 if session==1

gen session1_6=1 if session==6
replace session1_6=0 if session==1
	
*Figure 2 with ks test
/**Baseline***/

	   /*control vs. meal*/
ksmirnov pzbmi if child==1 & session==1, by(meal_control)  //p=0.928
local kspval_meal_control_1 = r(p) 
	di %6.3f `kspval_meal_control_1'

ranksum pzbmi if child==1 & session==1, by(meal_control)  //p=0.591
local ranksum_meal_control_1 = r(p) 
	di %6.3f `ranksum_meal_control_1'
	
   /*control vs. snack*/
ksmirnov pzbmi if child==1 & session==1, by(snack_control)  //p=0.368
local kspval_snack_control_1 = r(p) 
	di %6.3f `kspval_snack_control_1'
	
ranksum pzbmi if child==1 & session==1, by(snack_control)  //p=0.252
local ranksum_snack_control_1 = r(p) 
	di %6.3f `ranksum_snack_control_1'
	

twoway (kdensity pzbmi if child==1 & treatment==1 & session==1, scheme(s2mono) title(Baseline) ///
ylab(0(0.5)2)  legend(off) graphregion(color(white)) ytitle("") xtitle("BMI Percentile rank")) ///
(kdensity pzbmi if child==1 & treatment==2 & session==1) ///
(kdensity pzbmi if child==1 & treatment==3 & session==1)

	
/**After***/

	   /*control vs. meal*/
ksmirnov pzbmi if child==1 & session==3, by(meal_control)  //p=0.114
local kspval_meal_control_3 = r(p) 
	di %6.3f `kspval_meal_control_3'
	
ranksum pzbmi if child==1 & session==3, by(meal_control)  //p=0.045
local ranksum_meal_control_3 = r(p) 
	di %6.3f `ranksum_meal_control_3'
	
   /*control vs. snack*/
ksmirnov pzbmi if child==1 & session==3, by(snack_control)  //p=0.102
local kspval_snack_control_3 = r(p) 
	di %6.3f `kspval_snack_control_3'

ranksum pzbmi if child==1 & session==3, by(snack_control)  //p=0.013
	local ranksum_snack_control_3 = r(p) 
	di %6.3f `ranksum_snack_control_3'

twoway (kdensity pzbmi if child==1 & treatment==1 & session==3, scheme(s2mono) title(After) ///
ylab(0(0.5)2)  legend(off)  graphregion(color(white)) legend(size(small)) ytitle("") xtitle("BMI Percentile rank") )  ///
(kdensity pzbmi if child==1 & treatment==2 & session==3) ///
(kdensity pzbmi if child==1 & treatment==3 & session==3)

	
/**year 3***/
	   /*control vs. meal*/
ksmirnov pzbmi if child==1 & session==6, by(meal_control)  //p=0.215
local kspval_meal_control_6 = r(p) 
	di %6.3f `kspval_meal_control_6'
	
ranksum pzbmi if child==1 & session==6, by(meal_control)  //p=0.075
local ranksum_meal_control_6 = r(p) 
	di %6.3f `ranksum_meal_control_6'
	
/*control vs. snack*/
ksmirnov pzbmi if child==1 & session==6, by(snack_control)  //p=0.390
local kspval_snack_control_6 = r(p) 
	di %6.3f `kspval_snack_control_6'
	
ranksum pzbmi if child==1 & session==6, by(snack_control)  //p=0.095
local ranksum_snack_control_6 = r(p) 
	di %6.3f `ranksum_snack_control_6'
	

twoway (kdensity pzbmi if child==1 & treatment==1 & session==6, scheme(s2mono) title(Year 3) ///
ylab(0(0.5)2) graphregion(color(white)) legend(order(1 "Control" 2 "Meal" 3 "Snack")) legend(size(small)) ytitle("") xtitle("BMI Percentile rank") )  ///
(kdensity pzbmi if child==1 & treatment==2 & session==6) ///
(kdensity pzbmi if child==1 & treatment==3 & session==6)

